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ABSTRACT 



A classic problem in the study of the physics of sound 
in fluids is that of the finite amplitude plane wave propaga- 
ting in a viscous, unbounded medium. Though the solution to 
this problem is well known, for a given range of parameters , 
it would be desirable to develop techniques for solution of 
the problem which are not similarly limited. Through the ap- 
plication of the technique of parametric differentiation, 
this goal is realizable. This extension method transforms 
the governing non-linear differential equation to a linear 
equation in what may be termed parameter space; the equation 
is solved and a quadrature recovers the dependent variable 
solution. Parametric differentiation is conceptually straight- 
forward and has been applied in the past to a wide variety of 
non-linear equations. The method has b' n applied to the e- 
quation which describes finite amplitude plane wave propaga- 
tion in a viscous fluid in order to compare the results to the 
predictions of a known perturbation solution; the ultimate ob- 
jective being to utilize the technique for exact solution of 
the corresponding spherical and cylindrical wave problems. 

The first step in this process has been achieved; that is, 
numerical solutions for the plane wave case have been gener- 
ated for a range of values of the various viscous and non- 
linear parameters which are consistent with results obtained 
analytically. Further it is shown that solutions generated 
through the application of parametric differentiation may in 
fact have greater validity for certain ranges of viscous and 
non-linear parameters. 
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the coefficient of the first-order term in an assumed 
liquid pressure-density relation 



A~ , A~, A, - the variable coefficients of g , g , g , 

^ j t i X X L L X 

and g respectively 

XXL. 



arbitrary coefficients in a generalized boundary con- 
dition 



2 

small signal attenuation coefficient = vbk /2C for 
liquids, = vtH' + (y-1) /Pr] k^/2C Q for gases 



the coefficient of the second-order term in an assumed 
liquid pressure-density relation 



a viscosity number for liquids which represents shear 
viscosity and a phenomonological bulk viscosity 



the parameter of non-linearity = 1 + B/2A for liquids, 
= (y+1) /2 for gases 



C01, etc. - coefficients of the points of the finite-dif- 
ference cube 



small signal sound speed 
specific heat at constant pressure 
specific heat at constant volume 

interval or change in the variable which it precedes 
acoustic Mach number 



Neumann factor 



- 6 - 



n , n ' 


coefficients of shear and bulk viscosity respectively 


F 


the variable inhomogeneous term of the parameter space 
equation 


g 


parameter space variable = d^/dip 



9 V > 9 VV > etc. " notation for 9g/9x, 3 2 g/9x 2 , 9 3 g/9x 2 9t, 

X XX xxt XX 9 2 £/9x 2 , etc. 



9i- 5 i 


values of these variables for the ith value of ip, ip^ 


9 i.j 


indicial notation for the two-dimensional solution grid 


r 


an indicator of the strength of the non-linearity rela- 
tive to dissipation = 3ek/a 


Y 


ratio of specific heats = C^/C v 


i 

n 


nth order modified Bessel function of the first kind 


J 

n 


nth order Bessel function of the first kind 


k 


acoustic wave number 


K 


Bek product 


L [ ] 


a general linear operator 


X 


wavelength of acoustic signal, also used as an arbitrary 
real number in VonNeumann stability test 


m, n 


indicial notation as in g , corresponds to g . 

^m,n ^ irD 


N[ ] 


a general non-linear operator 



kinematic viscosity = n/p Q 

particle displacement amplitude = U /w 

particle displacement; subscript indicates the base 
solution, i=0 

an expression for a calculated particle displacement 
for output 

angular frequency 

constants in the fluid equation of state, for gases 
P = P Q , Q = 0; for liquids, P = P 0 C q 2 /Y/ Q = P - p Q , 

and y is determined empirically 
instantaneous pressure 
ambient pressure 
Prandtl number 

an arbitrary dependent variable, subscript indicates 
the base solution 

instantaneous density 

ambient density 

condensation = p-p /p 

o o 

a nondimensional distance = $ekx 



shock thickness 
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t , t 
' c 


temporal coordinate, characteristic time 


0 


an arbitrary complex number in VonNeumann stability 
test 


0 B. 

l 


boundary values of (p 


u, u 
' o 


particle velocity, subscript indicates base solution 


U 

o 


particle velocity amplitude 


X 


discontinuity distance = l/8ek 


X, x B 


a generalized spatial vector, the subscript indicates 
the boundary 


x ' x c 


spatial coordinate, characteristic length 


¥ 


viscosity number 


*' * o 


parameter of interest = vb/C X, subscript indicates 
base value 


V 


u/U c 


y 


a timelike variable in the Burgers equation solution 


c 


a non-physical variable used to transform the Burgers 
equation to a heat equation 
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I. INTRODUCTION 

The propagation of acoustic waves in fluids and certain 
crystalline substances has been the subject of long and ex- 
haustive study. Due to the inherent simplifications which a- 
rise in the examination of plane wave propagation, this prob- 
lem has received the closest scrutiny. In the general linear 
theory of sound propagation, the assumption must be made that 
particle velocity amplitudes are infinitesimal, allowing the 
reduction of the problem to the classic wave equation. How- 
ever, in many of today's acoustic systems this assumption is 
no longer valid, and solution of the fully non-linear wave 
propagation problem is required for analysis and prediction 
for practical systems. It has long been recognized that sound 
waves whose particle velocity amplitudes are not infinitesimal 
have phase velocities whose magnitudes change with the local 
particle velocity, a result which is not predicted by the lin- 
ear theory. This phenomenom received some attention from the 
philosophers of the nineteenth century, and with the exception 
of the efforts of two researchers of the pre-war decade was 
not treated until the last two decades at which time signifi- 
cant interest was generated because of the increasing sophis- 
tication and power levels of acoustic systems. 

The solutions of the 1930's form the basis of this re- 
search. Fubini [1] derived an implicit solution to the plane- 
wave problem in an inviscid fluid and reduced it to an explicit 
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solution for low Mach number. Fay [2] , the other researcher 
of the 30' s, solved the same equation with viscous effects in- 
cluded. Blacks tock [3] noted that though both solutions were 
entirely correct, they were restricted in their spatial re- 
gions of applicability. Applying weak-shock theory, he demon- 
strated the manner in which the two solutions are related and 
developed a single function which describes the solution for 
all space in the inviscid case. Keck and Beyer [4] performed 
a perturbation analysis which, for a specific range of para- 
meters, described propagation in the viscous case. Black- 
stock [5], using approximations which he demonstrated [6] to 
be equivalent to those used by Keck and Beyer, was able to re- 
duce the governing differential equation for plane wave propa- 
gation in a thermoviscous fluid to the Burgers equation for a 
boundary value problem, which yielded an analytic solution for 
all space and time. Propagation in relaxing fluids has re- 
ceived moderate attention [7,8], but has generally been a- 
voided. Blackstock [9] has derived analogous inviscid solu- 
tions for cylindrical and spherical waves, but as fate would 
have it, an analytic solution for cylindrical and spherical 
waves in a viscous fluid has not yet been obtained due to the 
fact that there exists no apparent simple transform, such as 
that utilized in the plane wave case, to obtain the Burgers 
equation . 

There are various parameters which describe the proper- 
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ties of the medium and two which relate to the source excita- 
tion, and taken together they give one an indication of the 
relative strengths of non-linear effects and viscous (dissipa- 
tive) effects. The parameters of the medium are v, the kine- 
matic viscosity; 4', the viscosity number; y, the ratio of spe- 
cific heats; Pr, the Prandtl number; 3, the parameter of non- 
linearity; and C , the small signal sound speed. The expression 

Y-l 

(4* + — ) is a collection of thermal and viscous coefficients 

which together with kinematic viscosity relate the thermovis- 
cous nature of gases. In liquids, such an expression is not 
readily determined and is generally replaced by a single term, 
b. 3 is equal to (y + l)/2 for perfect gases and to 1 + B/2A 
for fluids of arbitrary equation of state where A is the coef- 
ficient of the first order term and B/2 is the coefficient of 
the second order term in an assumed pressure-density relation. 
Analogous parameters for dielectric crystals allow them to be 
treated in a manner similar to that for thermoviscous fluids 
[10] . The two parameters related to the source excitation 
are U Q , the peak particle velocity, and to, the angular fre- 
quency of the source. Much of the previous effort in this 
field has made use of a parameter F = 3ek/a to express the 
relative strength of non-linearity and viscosity, where e = 

U o /C Q , the acoustic Mach number; k = w/C Q , the wave number; 

2 

and a = vbk /2C q , the small signal attenuation coefficient, 
r ~ 1 may be considered as the borderline for the inception of 
important non-linear effects [11] . 

Another important parameter is the discontinuity distance 
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X = ($ek) ■*"; at this point, in inviscid propagation, the slope 
of the particle velocity waveform becomes negatively infinite 
due to non-linear generation of harmonics, that is, X indi- 
cates the point at which a shock first forms. Shocks were 
first predicted by the inviscid theory, and the discontinuity 
distance which has become a reference point in most non-linear 
theory, was predicted by the Fubini solution, which, in fact, 
is not valid beyond the discontinuity distance. Shocks form 
because the finite amplitude of the particle velocity causes 
the compressive portion of the waveform to 'catch up' to the 
rarefactive portion of the waveform; because a multi-valued 
waveform is a physical impossibility, a shock must form unless 
dissipation reduces the wave amplitude sufficiently. Figures 
la and lb show graphically such a waveform. The utility of T 
has been extended by Fenlon [12] to indicate whether and where 
a shock will form in the viscous case; for plane waves T > 4.5 
indicates that shock formation will occur, that is, non-linear 
effects will be sufficiently strong, so that even with dissi- 
pation, shocks will form. Shock formation marks the end of 
the first of three identifiable zones of propagation of the 
finite amplitude wave; this is a region of strong non-linear 
effects. The second region of propagation is that in which 
non-linear effects initially dominate, but gradually taper off 
as absorption reduces the magnitude of the fundamental and 
even more rapidly the magnitudes of the non-linearly generated 
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harmonics. The limit of this region is reached when the rate 
of decay of the fundamental due to small signal attenuation is 
matched by the rate of decay due to the generation of harmon- 
ics, that is, when the particle velocity amplitudes are once 
again infinitesimal. 

The aforementioned solutions may be discussed with res- 
pect to their relations to the parameters of the problem. The 
Fubini solution predicts shock formation at the discontinuity 
distance, but the solution is no longer convergent beyond the 
discontinuity distance. Fay produced a solution for the most 
stable waveform, a sawtooth, to which the finite amplitude 
wave generates beyond about 3 discontinuity distances. The 
Fay solution contained r as a parameter, but did not reduce to 
Fubini for T °°; Blackstock's bridging function accomplished 
this. Keck and Beyer's perturbation solution is valid for 
very weak non-linear effects (relative to dissipative effects). 
Because of the limited number of terms generated, it is valid 
for T = 0(1) or less. Finally, Blackstock's Burgers equa- 
tion solution is valid for all r. Unfortunately the series 
which represents the steady state solution is very slowly con- 
vergent, limiting the usefulness of the solution. Finally, it 
should be noted that an inherent restriction on the size of 
the Mach number is assumed implicitly or explicitly in all 

I 

cases, with a maximum value of e = .1. Such a value would be 



stretching the applicability of certain approximations which 
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have been made in each of the above cases. 

In order to dispense with the aforementioned problems 
with the two viscous solutions, and specifically, to essenti- 
ally eliminate the approximations which were made therein, a- 
nother technique must be used which returns to the basic par- 
tial differential equation and solves it without approximation. 
Parametric differentiation is such a technique. Given a non- 
linear differential equation in which there exists a parameter, 
one need only assume that the solution is continuous in that 
parameter to a limiting value for which a known solution ex- 
ists. This technique has been applied in the past by Rappert 
and Landahl [13] to non-linear flow problems and by Harris [14] 
to aerodynamic sound produced by a rotating cylinder in a vis- 
cous medium. VanDyke [15] has linked parametric differentia- 
tion to work he has performed on the extension of perturbation 
series solutions. The application of parametric differentia- 
tion in this case to a problem whose approximate solution is 
known will (1) shed some additional light on the physics of 
the problem, and (2) demonstrate the applicability and utility 
of a new technique for use in studies of the propagation of 
finite amplitude cylindrical and spherical waves in a viscous 



medium. 
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II. THE EQUATION OF MOTION 

A plane wave propagating isentropically in a homogeneous 
viscous fluid may be completely described by the continuity 
equation, the Navier-Stokes equation for irrotational flow and 
an equation of state. The appropriate equations are: 



where p is the density, E, is the particle displacement, p is 
the pressure, n is the shear viscosity coefficient, g ' is the 
bulk , viscosity coefficient, and P and Q are constants dependent 
on the fluid considered. The subscripted parameters denote 
rest values of the respective quantities. For ideal gases, 

P = p Q , Q = 0, and y is the ratio of specific heats, C^/C^. 

For liquids, P = P 0 c 0 2 /Y' where C Q is the small signal sound 
speed, Q = P - p , and y must be determined empirically with 
the assumed equation of state above. The equation of state 
for liquids is often written in the form: 



1) Continuity p = (1 + 3£/8x)p Q 



( 2 . 1 ) 




3) State p = P(p/p o ) Y - Q [4] 



(2.3) 




(2.4) 



with = s the condensation. For small values of s, this 

Po 
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equivalent to the equation of state shown (2.3). This can be 
shown by writing 



p = P(1 + £ - £a ) Y - Q (2.5) 

Po 

and performing a binomial expansion for | p ~ p o | < 1, the result 

Po 

being 



P 



P + yP ( 



P-IP-O) 

Po 



y (y-1) 
2 



P ( P_£a) 2 + 

Po 



Q 



( 2 . 6 ) 



Then, with P - Q = p , A is equivalent to yP, and B to y(y-l)P 
and the ratio B/A = y-1. When one uses this form the implicit 
assumption is made that terms of °(( P ~ Pq ) an< ^ greater are not 
important. For any fluid of low compressibility this is 
clearly valid for small acoustic Mach number. (Note that p/p Q 
= (1 + If)' 1 and |f ~e) . 

Returning to Equations 2.1-3, p and p may be eliminated 
to yield a single equation in the particle displacement. From 
Equation 2.3 



9p _ Py ( _P_) Y 1 °!_P 
3x p 0 Po 9* 



(2.7) 



Substituting Equation 2.1 and its x derivative, 
3 r ~2 

(1 + T7- 2 -) , into Equation 2.7, one obtains, 



i£ = 

8x 



-P 



i!c 

o 3x 2 
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( 2 . 8 ) 



Substituting in the Equation 2.2, rearranging terms, and in- 
serting P = P 0 C 0 2 /Y/ one obtains the following equation: 



equation as it stands and apply parametric differentiation 
directly or make one additional simplification which resorts 
to the assumption that the displacements though finite are 
still not large enough to require the retention of cubic terms 
in particle displacement. Since previous effort has been dir- 
ected along the second path and since this previous work is 
the point for comparison to determine the validity of the ap- 
proach, this additional simplification is made. (1 + —■) ^ ^ 

dX 

is expanded binomially to yield: 



The first two terms only are retained since this factor multi- 
plies the right hand side. The resulting equation is: 




(2.9) 




(1 + » 







( 2 . 10 ) 




( 2 . 11 ) 
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where 8 = (y+l)/2 for gases and 8=1+ B/2A for liquids. This 
is in essence the equation treated by Keck and Beyer [4], a- 
long with the boundary condition £(0,t) = -5 0 cosiot, in their 
perturbation analysis, after terms of o((a/k) are eliminated. 
The Keck-Beyer solution is computationally simple, but it is 
severely limited by the labor required to calculate the terms 
of the perturbation series. In fact, according to the auth- 
ors, the perturbation series is probably only convergent for 
— 2 oix 

8ek(l-e )/4a < 1, which restriction corresponds, for exam- 
ple, to an acoustic pressure of .125 atmospheres at 1 MHz in 
water; one could easily argue that this pressure hardly qual- 
ifies as finite amplitude , but then if one is able to measure 
finely enough, any acoustic signal can be considered to be of 
finite amplitude. Actually, the series has been demonstrated 
to be everywhere convergent but is still limited in applica- 
bility, by truncation, to values of the parameters as defined 
above [7]. It is to be noted also that 8ek/4a = T/4 ; thus the 
Keck-Beyer solution will not allow T > 4, and since T > 4.5 is 
required for shock formation in the viscous plane wave case, 
this solution deals with non-shocked waves only. 

Blackstock has demonstrated that the governing equation 
may be reduced to a Burgers equation by applying the same 
assumption utilized in the perturbation analysis. This as- 
sumption, in fact, is the same one which allows one to arrive 
at an analytic steady state solution of the following problem: 
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3 2 £ vb3 3 £ _ ^2 3 2 £ 

’ 1 — — — — : — n — C . n~ 



3 1 2 3x3 2 1 o 3x 2 ' 



( 2 . 12 ) 



namely the small signal problem for plane wave propagation in 
a viscous fluid. This assumption, that a/k << 1, is valid for 
most fluids, with the exception of the most viscous ones, at 
any reasonable acoustic frequency. In any case, the effect of 
making this assumption is to allow one to utilize the expres- 



in reducing the full equation to a Burgers equation. This 
expression is, of course, the differential equation which des- 
cribes the propagation of a small amplitude wave in a lossless 
medium. The reduction procedure follows. Starting with Equa- 
tion 2.11 (an approximate one which has been obtained by as- 
suming a small Mach number), one first substitutes E, = “j^- 



sion : 



3 £ = _ _ 1 _ 9 £ 

3x C G 3t ' 



(2.13) 



o 



on the right hand side to obtain: 




(2.14) 



Then note the following identities from Equation 2.13: 




1 






(c) 



C 



o 
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xx 



— E, 

C o xt 



(e) 



’xx CTJ s tt 
o 



(f) (2 . 15a-f ) 



Substituting (d) into the viscous term of Equation 2.14, 



vb 



vb 



vb 



"2" Cvt- or> 2" ( or* 4 ^ttt *~0^xtt^ an( ^ 



C xxt 2C xxt xxt 2C 

o o o 



r ?+-+- + *» ( ^ttt C o ? xtt } 



’XX C z ^tt 2C 
o o 



r (2 - 16) 

o 



Substituting (e) and (f) into the non-linear term of Equation 
2.16, 



— r r = — (-F r -EE ) 
C H^xx C ' ^t^xx ^t^xx' 
o o 



$ / i r r _ i r t- \ 

C V C H^xt C 2 H*tt' 
o o o 



The equation after substituting and manipulating signs be- 
comes 



C 2 E 

oHt 



+ — (C 6 
2 ' o^xtt 



^ttt ) 



C 4 C 



= 6C (E E - c £J„J 



(2.17) 



Performing an operation which is the integral equivalent of 

3 3 

differentiating with respect to the operator -r— - C -r— , one 

dt O dX 

obtains 



C o?t - f Stt + C o«x - 



2 8 c o?t 



2 



(2.18) 
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Dif f erentiating with respect to ' t' and substituting u = 

C o u t " u tt + C o u x = 3C o u u t (2.19) 

Finally, with the variable transform 

x' = x t' = t - ^ 

o 

one obtains 

C o u x' " 3C o u u t' = I vbu t'f ' {2 - 20) 

which is Burgers equation for a boundary value problem. The 
significant point here is that one is considering an equation 
which is non-linear, and which includes viscous terms. The 
approximation a/k << 1, is applied two times to reduce the 
equation to Burgers form, each time to simplify viscous and 
non-linear terms. This is indeed an acceptable engineering 
approximation; however, it certainly is desirable to have the 
ability to solve the equation without resorting to this ap- 
proximation, particularly for highly viscous fluids and per- 
haps certain dielectric crystals with large stress coeffi- 
cients . 

As noted by Blackstock [6] , the Burgers equation solu- 
tion has one further limitation which might be considered more 
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severe. The steady state solution to the Burgers equation 
is 



Here, e is the Neumann factor (e =1, all other e = 2) : I 
n o n n 

is the nth order Bessel of imaginary argument; a = Bekx, a 
non-dimensional space variable; and y = wt', a non-dimensional 
temporal variable. For values of r > 50, the series represen- 
tation of the solution is very slowly convergent, and there- 
fore difficult to use practically in numerical analysis. Con- 
sequently, for T > 50, Blackstock utilized an asymptotic ap- 
proach to the solution, which in fact is not valid near the 
origin. To solve the Burgers equation for Bekx << T, he re- 
sorts to a perturbation analysis in T \ but recalling T= Bek/a, 
this is in some sense similar to the Keck-Beyer approach which 
utilized e as the perturbation quantity to solve the same basic 
equation. Blackstock' s solution for small x contains terms of 
first order in the perturbation parameter; for r >> c, this 
should indeed be appropriate, however the coefficients derived 
contain an infinite series of Bessel functions which may be 




cos ny 



( 2 . 21 ) 



where £ is related to u by, 



u = U o f- [log C] y 



( 2 . 22 ) 
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slowly convergent themselves. Herein lies another motivation 
for the parametric differentiation solution, the desire to 
produce a unified solution applying the same assumptions and 
having the same limitations throughout. 

In any event, there is no simple transform which will re- 
duce either (1) the full equation valid for high Mach number, 
or (2) cylindrical and spherical finite amplitude wave equa- 
tions to the Burgers form. Therefore, an alternative ap- 
proach to the problem is required. 
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III. PARAMETRIC DIFFERENTIATION 

Parametric differentiation is a procedure which allows a 
non-linear equation which involves a parameter to be solved by 
transforming the equation to a linear equation in parameter 
space. This is particularly valuable for cases in which nu- 
merical solution will be required because it eliminates the 
possibility of obtaining multiple roots in solving non-linear 
difference equations. It is essential to note at this point 
that some insight into the nature of the solution at the ex- 
tended range of the parameter of interest would be extremely 
helpful. If the solution is singular in the parameter, that 
is, if the nature of the solution changes radically as one ap- 
proaches the limiting value of the parameter in the known base 
solution, then the technique is not applicable. Thus, para- 
metric differentiation would never be applicable to the entire 
class of singular perturbation problems. The method of appli- 
cation follows. 

Suppose one is given the general problem: 



N[(j)CXf t; ip ) ] = 0 



(3.1) 



<|>[(X B , t; ip)] 



0 B (t; i|») 
1 



(3.1a) 




(3.1b) 



acp ( (X B r t; ip) ] + c<j, n (x B , t; ip) 




(3.1c) 
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where N is a non-linear operator; X is a spacial vector coor- 
dinate; t is a time coordinate; if; is the parameter of inter- 
est; and the 0 r , 's are functional values of 4, d> , or a6 + cd) 
at X D as appropriate to the problem, n indicating the normal 
derivative. Initial values could also be prescribed if appro- 
priate. If a solution <p Q = <f>(X,t; exists at some limiting 

value of the parameter if; Q , and it is desired to obtain a solu- 
tion to the same problem with identical boundary and/or ini- 
tial values at + Aip, then parametric differentiation may be 
applied. First, differentiate the governing equation and 
boundary and/or initial conditions with respect to the para- 
meter to yield: 



L [g (X, t; ip) ] = 0 



(3.2) 



g(}f B , t; if; ) = g B (t; if;) 



(3.2a) 




(3.2b) 



ag (Xg , t; if;) + cg n (X fi , t; if;) 



g B (t; if;) 



(3.2c) 



where 




(3.3) 



One should realize that this procedure cannot be applied if 
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the parameter appears in the basic equation raised to a power 
less than one, and if the limiting value of the parameter \p Q 
= 0 . 

The operator, L, is a new linear operator obtained through 
differentiation. The new equation with associated boundary 
and/or initial conditions may now be solved in parameter space, 
either analytically or numerically, for successive values of 
the parameter ip . As the equation is solved for each value ip , 
a quadrature is performed which solves Equation 3.3 and re- 
turns a value or functional form which, when applied to the 
solution for yields the desired solution at iJk . There 

is no numerical restriction on the range of ip; however, there 
certainly may be a physical limit, and some insight is re- 
quired here to avoid blind application of the method to obtain, 
for instance, solutions for a range of ip which is beyond the 
region of validity of the basic equation. The quadrature may 
be expressed as follows: 

_ fi 

< P i (X, t; ijj i ) = (X , t; + I gcjijj (3.4) 

♦i-1 

For more complicated equations which require numerical solu- 
tion, this step can be quite complex, depending on the form of 
the g versus \p curve. In some cases, trapezoidal integration 
may be sufficient; in others, higher order integration techni- 
ques may be required. It is clear in examining the form of 
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the equation above, that this is, in essence, a recursion type 
relation and consequently it demands a known initializing 
function for i=l. 

In this work, the initializing function is the Fubini 
solution to the problem, 



(1 + 



3 £, 26 

3x ; 
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= C 



o 3x 2 



(3.5) 



with the boundary condition £(0, t) = - E o cosu)t, this being a 
duplicate of Equation 2.9, without the dissipative term. Thus, 
the parameter of interest here will be some non-dimensional 
quantity which contains v, and the quadrature will be per- 
formed with respect to the non-dimensional quantity. An ex- 
plicit form of the Fubini solution was given by Keck and 
Beyer [4] in terms of particle velocity, 

<» J (ntcx) 

u ( x , t) = 211 l n „ - — sinn (tot - kx) (3.6) 

n=l 



where J is the nth order Bessel function and k = Bek. This 
n 

solution is limited in Mach number (^ << 1) since a binomial 

L o 

expansion in Mach number was utilized to arrive at this result. 
An exact implicit solution is available which could be used for 
analysis at greater Mach numbers, 
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However, the numerical problems associated with an iterative 
solution at many points would be difficult to resolve. This 
solution predicts the formation of the shock at the point 
where the velocity profile becomes negatively infinite, that 
is, where 3u/3x -> -°°. Differentiating Equation 3.7 with res- 
pect to x. 
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This quantity becomes infinite when 

_2y_ 

x = (1 + Iz1^)V- 1 /1+1Uq « cos( , ,3.10) 

o o o 

Noting that 8 = — e = ■— = k, and that the first shock 

o o 

will occur at the head of the waveform where cos ( ) =1, and 
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also that u = 0 at the shock, the discontinuity distance X 
= -rij- is obtained. X is the limit of convergence of the ser- 

pGK 

ies solution stated above. Thus, using the Fubini explicit 
solution as a starting solution, one is limited to application 
of the technique of parametric differentiation in the region 
0 < x < X and for Mach numbers which are small compared to 
unity. It is clear from this discussion that one is fundamen- 
tally limited when applying parametric differentiation by the 
quality of the base solution. It is unfortunate that in this 
case, the base solution is limited with respect to Mach number; 
on the other hand, signals of sufficiently large amplitude may 
no longer be thought of as acoustic waves and may propagate in 



a different manner. 
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IV. APPLICATION OF THE METHOD 

It is one thing to discuss the use of parametric differ- 
entiation in abstract terms, and quite another to execute its 
application, there being no real 'theory' with respect to how 
it should be applied. In any case, the first step in the sol- 
ution technique is similar to that for most other numerical 
solution procedures, that is, non-dimensionalize the equation. 
Starting with the reduced equation 



€ 



xx 



1 p , Vb r 

c 1 Ht c z Scxt 

o o 



2$r £ 

X XX 



(4.1) 



one introduces characteristic lengths and times. This being a 

wave propagation problem, it is reasonable to select x c = A and 

t = ■— , where A is the wavelength; then the non-dimensional 
c o 

variables are; 



x = x/x^ = x/A t = t/t c = tC Q /A and E, = K/* c = eJx 



Then, Equation 4.1, in terms of dimensionless variables, is: 



A /s ^ . y_ 

+ .r. 

s xx ^tt C A 

o 



xx t x xx 



(4.2) 



This equation is as one would expect for a 'good' non-dimen- 
sionalization ; the coefficients of the first two terms, which 
in fact should dominate, are of 0(1). Alternatively, with the 
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equation in the form 



(1 + «*> 26(e tt - 



) = C?£ 



O XX 



(4.3) 



a similar non-dimensionalization will yield 
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which demonstrates that the relatively large numerical value 
of the coefficient of the non-linear terms in 4.2 is simply a re- 
sult of the expansion performed, and should not be considered im- 

vb 

proper. The viscosity coefficient appears as p — r- which is 

c o A 

like an inverse Reynold's number. However, because C Q is a 
propagation speed and not a particle velocity, this is not a 
classic acoustic Reynold's number. In any case, this quantity, 
which will be termed ' \p ' , is the parameter of interest, on 
which the quadrature will be performed. The solution for i^=0 
is known; it is the Fubini solution. Then one has 



/ \ 

u (X, t; ) = 2U T ^n (nKx) sinn (cut - kx) (4.5a) 

o r o o L , n«x 

n=l 

or, in terms of particle displacement 



£ (X, t; 4* ) = -2^ y n ( 2 nKx ) cosn (cot - kx) (4.5o) 

^o r o w n kx 

n=l 

which corresponds to the function <j> o in the theoretical devel- 
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opment of parametric differentiation, and with H = —2. 

O Wo 

The next step in the development is to differentiate E- 
quation 4.2 with respect to ty , yielding 
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" g tt + 






xxt 



£ = 2BC g 

xxt x^xx 



+ 



2 $£ g 

XX 5 X 



(4.6) 



where 



g 




(4.7) 



and the umlaut notation has been dropped from the non-dimen- 
sional variables. One would immediately look at this equation 
and express horror at what has been done; there seems to be 
added complexity rather than simplification. However, one can 
now treat the equation as an inhomogeneous, linear equation 
with variable coefficients, which should make it more tract- 
able for numerical solution. The only non-linearity which re- 
mains is found in Equation 4.7. It should be noted that the 
non-dimensionalization in terms of the size of the coefficients 
has not been adversely affected. In fact, the dominance of the 

first two terms has been enhanced since £ , £ ~e and e << 1. 

x xx 

Alternatively, differentiating Equation 4.4 one obtains 



(2BH1 + 5 x ) 2B - 1 g x (C tt - « xxt > ♦ (1 - 5 x ) 28 (g tt 



^xxt ^xxt^ g xx 



(4.8) 
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which really looks as if this "simplification" could have been 
done without. However, if one multiplies through by (1 + E, ) , 

X 

the following equation is obtained 

(26) (1 + 5 x ) 26 (c tt - *S xxt >9 x + (1 + C x ) 26+1 (g tt - *9 xxt 



- c .) = (l + c )g 

xxt ^x ^xx 



or 



+ (1 + ? x ) 2B+1 (g tt - *g xxt - 5 xxt ) = d + C J g 



X 



X J xx 
(4.9) 



which is of the same form as Equation 4.6. This now is the 
equation which will be utilized to solve the finite amplitude 
plane wave propagation problem without approximation, save 
those which may be necessary to obtain a base solution. In 
this case, this allows a significant extension because this 
equation does not require the approximation a/k << 1, and 
therefore, the Fubini solution may be used as a base solution 
for the problem of propagation in highly viscous fluids since 
it is limited by Mach number only. 

The quadrature required to extract the desired solution 
is essentially dependent on the problem itself. Landahl and 
Ruppert [13] employed Runge-Kutta integration to solve a boun- 
dary-layer problem. In many cases, it is not necessary to re- 
sort to such exotic techniques; in this analysis, g varies quite 
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slowly for small values of \p; a rapid non-linear change in g 
with \p only occurs as ip becomes very large, that is, as the 



viscous ones, a value of ^ sufficient to produce rapid varia- 
tion of g with ip would not be attained at any reasonable acou- 
stic frequency. Consequently, the integration technique em- 
ployed in this analysis was a very simple straight line inte- 
gration, which required one iterative step. The initial sol- 
ution of the equation for g on the ith step was computed based 
on the ith value of ip and a solution £^(x, t) which was ob- 
tained in the following manner 



With this new solution, g^(x, t) , a modification of £^( x , t) 
was performed as follows: 



viscosity becomes high. For most fluids, including very 



(x, t) 
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(x, t) = £. (x, t) + i(g. (x, t) 

1 2 1 1 






( 4 . 11 ) 



The net result being: 




( *i - h-i’ 



( 4 . 12 ) 
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which is in essence a trapezoidal integration. Clearly, this 
scheme would not be useful for g functions which vary rapidly 
and non-linearly with ip , that is, if Equation 4.7 exhibits 
highly non-linear behavior. This scheme could be used in an 
iterative manner, however, with continuing refinement of the g 
values at each step. The entire process used to arrive at the 
extended solution for £ (x , t) is shown schematically in Figure 
2. It may be noted that there is a requirement to obtain var- 
iable coefficients at each step, this of course is a situation 
which is not enviable in a numerical solution in which the co- 
efficients cannot be described analytically, particularly when 
these variable coefficients involve derivatives of the depen- 
dent variable. In this case, all the coefficients involve de- 
rivatives. A simple solution to this problem has not been 
found; however, in this analysis, numerical differentiation of 
the E, solution to obtain the coefficients does not seem to 
have affected the results. 

Another not so simple problem in the application of para- 
metric differentiation is the treatment of the boundary condi- 
tions. One must be extremely careful in not trying to read 
too much into the problem. In this analysis, two approaches 
were considered for determination of boundary values for appli- 
cation in the finite difference scheme, the result being that 
the most natural and straightforward approach was eventually 
chosen as the proper one. In analyzing this problem, one can 
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make some subjective evaluations of the expected behavior of 
g based on what the gAi|> product represents. Dissipation will 
tend to slow the non-linear generation of harmonics as the 
wave propagates, and if dissipation is strong enough, one 
would expect the waveform to remain sinusoidal, if such is the 
nature of source excitation, no matter how strong the non-lin- 
earity is in absolute terms, that is, for F < 1 the waveform 
should not steepen appreciably. In any case, near the bound- 
ary the effect of gA\p should be to decrease the amplitude of 
the signal. A\p is always positive, therefore to obtain such 
a reduction in amplitude, one would expect the sign of g to be 
opposite that of the base solution. Attempting to produce an 
analytical model which fulfills one's expectations near the 
boundary may be self defeating. In this analysis it was ini- 
tially considered desirable to introduce a small perturbation 
of the boundary condition in terms of \p , and consequently the 
boundary condition was rewritten as, 

£(0, t) = - j cos(cot) (4.13) 

which yielded the intuitively desirable result in terms of 
sign. 



3 



(Uo, 



t) 



g(°/ t) = — cos (cot) 



(4.14) 



This model also seemed appropriate since it produced a bound- 
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ary value of g which decreased slightly as i/j increased, corres- 
ponding to an anticipated increase in the viscous effects. If 
one recalls the nature of the quadrature, the only effect of 
this model would be to force the solution near the boundary to 
take on the desired form. This 'forcing' of the solution by 
injecting one's subjective expectations into the modelling of 
the boundary conditions was not particularly bad, but neither 
did it accomplish anything with respect to improving the solu- 
tion. If one instead blindly applies the — operator to the 



aijj 



boundary condition, the result is 



g (0 , t) = t^-(-E o cos (wt) ) = 0 



( 4 . 15 ) 



which does absolutely nothing to convince one that the proper 
sort of quadrature will occur near the boundary, but instead re- 
lies on the differential equation to properly predict the g 
values everywhere. In the final analysis, this was the cor- 
rect procedure; there were in fact very small differences in 
the solution values obtained with the forced boundary condi- 
tions, and those obtained with the 'natural' condition. 

One may wonder why initial conditions are not important 
in this problem since it is essentially hyperbolic in char- 
acter, and such problems are generally thought of as initial 
value problems. The answer is that when one is interested in 
the steady state solution only, for a problem with a periodi- 
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cally varying boundary condition, and in this case that is the 
solution desired for £ (x , t) , one may set the condition that 
the initial values occurred at some sufficiently large time in 
the past that all start-up transients have disappeared, and 
therefore only the boundary conditions determine the behavior 
of the solution [16]. It is only natural to expect this physical 
interpretation to be carried over into the parameter space ap- 
plication, and consequently, initial values have been of no 
concern. 

It is valuable to note that the final result for the cho- 
sen maximum value of ip need not be the only output. For a 
given Be product, a solution for any T consistent with the ac- 
curacy obtainable (as T ■+ <», the accuracy required to obtain 
meaningful results is increased for a given Be because the 

differences between the inviscid and viscous solutions go to 

6 £ 

zero) , may be obtained by selecting the appropriate xp = ^ and 
simply calling for output at the end of the second iterative 
step. For a given fluid this would allow one to analyze the 
losses at different frequencies for a given Be and, in this 
case, out to the discontinuity distance. In fact, tabulations 
of g versus ip at specific points could be stored so that the 
solution could be extracted in the future simply by integra- 
ting and applying the correction at the appropriate point in 
the base solution. 

It is useful to compare the base solution waveform with 
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the anticipated extended solution waveform, the comparison in 
this case being very simple because the extended solution or 
a very good approximation thereto is well known. As it de- 
velops, the incremental differences between the base solution 
and the extended solution are expected to be in phase with one 
another, if they are written in the form 

A £(x, t) = £ f (x, t) - £ q (x, t) (4.16) 

where A£(x, t) is the difference between the calculated solu- 
tion and the base solution when integrating Equation 4.7 from 
ij ) q to The observance of A£(x, t) varying in the expected 

manner is a key clue to the success of the application of the 
method. This behavior was in fact observed in this analysis, 
and will be discussed in a following section. Once one has 
developed an equation in parameter space and made a subjective 
evaluation of what he expects to be the functional form of g 
and of A E, (x, t) , he is ready to proceed to the numerical tech- 
nique to be used for solution. Ideally, it would be hoped 
that the transform to parameter space would lead to an anal- 
ytic solution of Equations 4.6 or 4.9 and Equation 4.7, or at 
the very least, Equation 4.6 or 4.9. Such is not the case 
here, and numerical solution of both was required. Thus, to 
analyze the results, and prove the worth of the method, a sub- 
jective feel for the results was required. 
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V. _ _ FIN ITE DIFFERENCE EQUAT IONS 

After the equation of motion has been derived, non-dimen- 
sionalized and its analog in parameter space has been obtained 
through differentiation, it may still be a significant task to 
obtain the solution of the linear g equation. In this case. 
Equation 4.6 or 4.9 describes the function g with boundary 
condition g(0, t) = 0. In analyzing the equations, it is noted 
that they are essentially wave-like in nature, at least close 
to the boundary, since the g and g terms dominate. Adopt- 

XX L. L. 

ing the approximation E, ~ E,. (in non-dimensional form) and 

X u 

differentiating with respect to x, one notes that E, » -u 

XX X 

and since u -*•-«> at the shock in the inviscid case, then E, 

x xx 

-»■ 00 at the shock. This wave-like behavior may not then be so 
readily observable away from the origin since the coefficients 
involving E, may become very large in the regions in which 

• XX 

the shock is forming, or has formed in an extended solution, 
in which case the g term becomes important. The g . term 

X XX L. 

remains small since ip is generally much less than one. The g 

X 

term, of course, is the one which results from the inclusion 
of non-linear terms, and it is entirely consistent that it be- 
comes important in the regions of the waveform which exhibit 
the result of non-linear propagation, namely in the region at 
the head of each compressive portion of the waveform and at 
the tail of the rarefactive portion. Once a shock has formed, 
it is conceivable that the g term will alter the character of 
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the equation to the extent that it is no longer wavelike; this 
is not of immediate concern since this analysis is restricted 
to the first region of propagation. Therefore, in attempting 
to solve this equation, it was treated as if it were wavelike 
throughout. As previously noted, the literature is rich with 
various methods of solution of the fundamental equation of mo- 
tion, techniques which, if applied judiciously and with essen- 
tially the same approximations as those utilized to solve the 
equation of motion, could possibly be applied to Equation 4.6 
or 4.9 to yield analytic solutions. However, because of the 
form of the variable coefficients, this analytic solution 
would at best be extremely difficult and most likely would be 
impossible; so in this analysis solution, by finite differences, 
of the parameter space equation was accomplished, and numeri- 
cal integration of Equation 4.7 was performed. Since the e- 
quation is essentially hyperbolic in nature, a finite differ- 
ence scheme appropriate to a hyperbolic system was the obvious 
choice. 

There are two broad categories of finite difference 
schemes for partial differential equations, which are termed 
explicit and implicit. Adopting a standard notation for a 
nine point difference mesh, Figure 3, an explicit scheme, is 
one which has only one unknown in the i+1 column or the j+1 
row, depending on whether one is "marching" in the i or j dir- 
ection. An implicit scheme is one for which there are two or 
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more unknowns in the i+1 column or j+1 row. To solve the ex- 
plicit scheme for a single value in the future (marching dir- 
ection) one need know only the values of the dependent varia- 
bles at all present and past positions on the mesh. The terms 
future, present and past apply in a spatial as well as tempor- 
al sense. Thus, for the explicit scheme, if one has initial 
conditions valid for all space, or if one has initial condi- 
tions for a half space and a boundary condition, a complete 
solution can theoretically be generated numerically for a pro- 
perly-posed problem. The implicit scheme requires solution at 
two or more points based on present and past points. Thus, 
one generates two unknowns in the first application of the 
mesh, and an additional unknown with each subsequent applica- 
tion, yielding N+l unknowns with N equations. Thus, an addi- 
tional condition is required in order to apply an implicit 
scheme; generally it is a boundary condition at some point in 
space. In any case, this requirement essentially restricts 
the use of implicit schemes in infinite domain problems such 
as this one. 

Having been limited to the use of explicit solution 
schemes, there exists one additional problem of significant 
import; that is, stability of the scheme. Unfortunately, it 
is difficult to consider stability until one has actually de- 
veloped the finite difference equation, stability of the scheme 
being a function of the way the scheme itself is developed. 
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Proceeding to finite difference approximations of the Equation 
4.6 or 4.9, one finds a literature which is rich in references 
to the wave equation for infinitesimal amplitude waves, in a 
lossless medium, but which has very little in the way of ref- 
erence to non-linear equations or equations for a viscous 
medium. Recalling that the parameter space equation is essen- 
tially wavelike as long as the coefficient of the g term is 
relatively small, one can use the finite difference approxima- 
tions to the wave equation as a starting point. Since an es- 
sential aim of this study was to demonstrate the utility of 
the method of parametric differentiation in non-linear wave 
problems, a second order finite difference scheme was selected; 
it is clear that greater accuracy could have been obtained with 
a higher order scheme. The following finite difference approx- 
imations were used for each of the derivatives of g in either 
Equation 4.6 or 4.9: 
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where Ax and At denote the step sizes in space and time res- 
pectively. It must be stressed that these approximations, par- 
ticularly Equations 5.3 and 5.4, are not the only ones which 
could be used. For instance, a forward or backward difference 
could have been used in place of the central difference of E- 
quation 5.3. This would, however, have changed the order of 
the approximations. The form of Equation 5.4 is that of a 
second central difference in space and a backward difference 
in time, the backward difference being necessary in order to 
preserve the explicit nature of the scheme. A central differ- 
ence or forward difference in time for Equation 5.4 would have 
not only destroyed the explicit nature of the scheme but would 
also have yielded an unknown to be determined at the point 
g^ + ^ whose coefficient, ip , would have been much smaller 

than those of the other terms, and could have led to an in- 
stability, even if one were using an implicit scheme. This 
problem could be resolved in an implicit scheme by substitu- 
ting for Equation 5.1 a weighted time average of second order 
approximations of g 

X X 

Representing the variable coefficients of g r , g , , g , 

XX U u X 

and 9 xxt by A^, , A^, and A^ respectively. Equation 4.6 or 

4.9 becomes, upon substituting the finite difference approxi- 
mations 5.1-4: 



~ A 4*^i-l j_]_/(Ax) 2 At + (A-^/ (Ax) 2 + A 3 /2 Ax + A 4 /(Ax) 2 At) 
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* g i-l j + ( -A 2 /^ At ) 2 + 2A 4 / (Ax) 2 At) ^ ^ + (-2A ]L /(Ax) 2 
+ 2A 2 /(At) 2 - 2A 4 / (Ax) 2 At) •g i ^ - A 2 * f j +1 / ( At ) 2 
- A 4'9i + i j_ 1 /(Ax) 2 At + ^/(Ax) 2 - A 3 /2Ax + A 4 /(Ax) 2 At) 



i+1 , 



F 



(5.5) 



where F is the inhomogeneous term. At this point, one may 
look at Equation 5.5 and note that it has one unknown in the 
j+1 row if one is marching in space, and two in the i+1 column 
if marching in time; so that it is well suited to treatment of 
an initial value problem explicitly, but is implicit if used 
for a boundary value problem. In order to apply the scheme in 
an explicit manner to the boundary value problem, one must make 
an additional approximation. To solve for g(i+l,j), one needs 
information in the i-1 and i columns, plus information at the 
point g^ + ^ j-l' Because the third order derivative is present, 
there is no way to avoid this problem. However, considering 
the fact that the coefficient of the third order derivative is 
small compared to the coefficients of the other terms, it would 
not be unreasonable to presume that if one had an approximate 
value of g^ + ^ j-l' ^hen a solution of reasonable accuracy could 
be obtained for g^ + ^ j, which would be based on five pieces of 
accurate information and a single piece of approximate infor- 
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mation. Considering Figure 4, which depicts the finite-dif- 
ference grid with the coefficients of the various points in 
the solution cube annotated thereon, one notes that the ap- 
proximation required at g^ + -^ need be made only once in 

each spatial step. After solving for g^ + ^ j at the 'base' of 
a single spatial column, g^ + ^ is now known for all subse- 

quent applications of the cube as one marches out in time. The 
approximation used to arrive at this "corner point" is a Tay- 
lor's series expansion in which the derivatives at the point 
about which the expansion is taken are approximated by back- 
ward differences accurate to O(Ax). Again greater accuracy is 
obtainable, but was not required because of the relative small- 
ness of the coefficient. The resulting expression for the 
"corner point" is, 



i+1 t j - l 



= 2. 5*g 



l » 



j-l 



- 2 



* g i-l,j-l + ^i_2 , j -l 



(5.6) 



It must be noted that this approximation is a potentially weak 
link in the method, due to the fact that near the head of the 
waveform the x derivatives become large, and a truncated Tay- 
lor's series may not be of sufficient accuracy. In any case, 
it did not appear to significantly affect the nature of the 
results in this analysis. 

In the discussion of explicit finite difference schemes 
one must inevitably be concerned with the stability of the 
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scheme as it relates to step size. According to the Courant- 
Friedrichs-Lewy criterion for stability of finite difference 
schemes [17] , the ratio of the step sizes must be such that 
one does not attempt to find the solution at a point which is 
outside the region of influence of the known data, the region 
of influence being determined by the characteristics emanating 
from the endpoints of the known data. Unfortunately, this 
problem is difficult to address here because of the nature of 
the governing differential equation. Typically, third-order 
partial differential equations are not treated as such, but 
are instead reduced to some more tractable form; in this case 
the approach using the Burgers equation was such a reduction. 
Consequently, little appears in the literature with reference 
to stability for third-order equations. An unsuccessful in- 
vestigation was conducted to attempt to arrive at the exact 
characteristics by reduction to canonical form. However, this 
equation is wave-like, and one could reasonably expect the 
characteristics to be similar to those of the wave equation, 
and that the stability criterion would be similar in form. 
Working with this foreknowledge and applying a modified Von 
Neumann [17] stability test, an exact stability criterion will 
be established. The VonNeumann test essentially consists of 
examining possible exponential solutions to a finite differ- 
ence scheme to determine when they may grow without bound 
given finite initial or boundary conditions. With this know- 
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ledge one may adjust the stepsize so that the scheme will be 
stable. Unfortunately, the VonNeumann test is applicable to 
finite difference equations with constant coefficients; in 
this case the coefficients are variable, and the only recourse 
is to treat the coefficients as if they were constant, derive 
the stability criterion, and then discuss the implications of 
their variability with respect to how it affects the stability 
criterion. 

The exponential solutions to the scheme which may grow 
without bound will have the form, 



g 



m, n 



im0 

e 



inA 

e 



(5.7) 



where the subscript notation m,n corresponds to the notation 
i,j used previously. A is any real number and 9 is an arbi- 
trary complex number. Thus, boundary data may be written, 
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inA 

e 



(5.8) 



and the right hand side may be considered as a typical term in 
a Fourier expansion of the boundary data. Substituting Equa- 
tion 5.7 into 5.5, ignoring the inhomogeneous term since it 
will not affect stability, and dividing through by e im ®e ln \ 
one obtains 

-A^e ^e (Ax) 2 At + (A^/(Ax) 2 + A^/2Ax + A 4 /(Ax) 2 At)e ^ 
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+ (-A 2 /(At) 2 + 2A^/ (Ax) 2 At) e ^ + (-2A^/(A x) 2 + 2A^/ {^) 



- 2A^/ (Ax) 2 At) - / ( At) 2 - A^e^e (Ax) 2 At 



+ (A.. / (Ax) 2 - A /2 Ax + A./ (Ax) 2 At) e 10 = 0 (5.9) 



Rearranging terras such that they are grouped by coefficients, 
the following results. 
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. ,e +e -2. ,e +e -2. 

A 1 ( (Ax) 2 ) A 2 ( TAtp ) 



i0 -i0 

tv ( S : — ) 

3 l 2 Ax ' 



. . ,e i0 + e- ie -2 w l-e“ i \ „ 
A 4 ( (A xP } ( ' At } ” 0 



(5.10) 



Recognizing the various terms as trigonometric forms, one ob- 
tains 



(4A^sin 2 -^-) / ( Ax) 2 + (4A2sin 2 -^-) / ( At) 2 - iA^sinO/Ax 



A^ (4sin 2 |-) (1-e ^ ) / ( Ax) 2 At = 0 



(5.11) 



Now solving for the imaginary portion. 



-A^sinO/Ax - 4A^sin 2 |- sinX/(Ax) 2 At = 0 



which identity indicates that the product, 



(5.12) 
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(5.13) 



should hold, but this is clearly absurd since it indicates 
that as Ax -* O, At may become arbitrarily large as well as 
negative. Reason and foreknowledge of the stability criterion 
for the wave equation dictate that this result has no meaning. 
Returning to Equation 5.11 and equating the real parts, one 
obtains 



Now one argues that to prevent the assumed exponential solu- 
tion from increasing without bound, 0 must have no negative 
imaginary part. However, because of the quadratic form, if 0 
has imaginary roots, they will appear in conjugate pairs; 
therefore 0 must be real. Since A was assumed arbitrary, this 
requires the following restriction, 



(-A^sin 2 -^-) / (Ax) 2 + (A^sin 2 ^-) /At 2 - (A^sin 2 -|) (1-cosA) / ( Ax) 2 At 



0 



(5.14) 



or 



• 2 9 

ln 2 “ A x (At) z +A 4 At (1-cosA) 




(5.15) 
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(5.16) 
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Now, noting that the maximum value of the expression above is 
obtained when 1-cosA = 0, it is clear that a sufficient ex- 
pression for the above is, 

A, (Ax ) 2 

0 i A^iAty^ - 1 (5 - 17) 

which is of precisely the same form as the stability criterion 
for the wave equation, in which case A^ = = 1 . In this 

2 8+1 

case A, and A„ are variable with A, = 1 + E, and A 0 = (1 + £ ) p 

12 1 v x 2 v x 

for Equation 4.9. Since E ? is approximately equal to the lo- 

X 

cal acoustic Mach number, which is taken to have a value of 

-2 

0(10 ) in this analysis, the criterion is only slightly dif- 

ferent from that of the wave equation. However, it is impor- 
tant to realize that the ratio of step sizes cannot be unity 
as one may have anticipated. It is also comforting to note 
that A_, the coefficient of the g term, did not enter into 

J X 

the criterion in any way, lower order terms generally having 
little effect on stability. This statement, however, deserves 
further comment. Recalling the discussion of the size of the 
respective terms in the parameter space equation, the g term 
has a coefficient, £ , which will become very large in the 

region of a shock. It is inconceivable that this term will 
not affect the stability of the scheme in such a case. The 
point is that when £ becomes very large the entire nature of 
the equation does change with the g term being the critical 

X 
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term. In this instance the stability analysis is no longer 
valid, and one would be forced to rely on some other stability 
test. This problem is not great here because the analysis is 
restricted to the first zone of propagation. Also, the dis- 
cretization of the displacement data limits the maximum value 
that £ can attain in the numerical scheme. This is, as yet, 
an unresolved numerical problem; to allow £ _ to attain its 
true value in a nearly inviscid fluid, an exceedingly small 
step size in space would be required. However, if one consid- 
ers the actual size £ may attain in a typical propagating 
wave, this problem may be put in a proper light. Blackstock 
[5] has defined a shock thickness for the viscous problem, 

T = 2Acosh ^/tt/A (5.18) 

where T is the thickness, defined as the spatial distance from 

trough to peak of the waveform, and A = 2(1 + 0 ) /tiT . As an 

example, consider the case when T = 100 and a = 1; then T = 

8 *“ 1 

TnrT - cosh 5 it ^ .08. In this case, £ changes from a maximum 
IOOtt ’xx 3 

negative value to a maximum positive value over a change in x 
of .08. Since £ -u , the change in u over this distance 

XX X 

is approximately twice the Mach number, which as noted before 

_2 

is taken to be of 0(10 ). Then an approximate value of £ xx 

may be taken as Au/Ax ^ .25; therefore, while the g term does 
become important, it does not dominate. For larger values of 
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r, it will become increasingly important, will in fact domi- 
nate the equation, and will require smaller step size near a 

shock to correctly estimate £ 

J xx 

Another consideration in the finite difference scheme is 
the handling of boundary values. Treated as an initial-value 
problem, there would be no problem in starting the solution, 
since £ and would be everywhere prescribed for t = 0. In 
this case, only one boundary condition is available for x = 0, 
the other condition being that £ -*■ 0 as x -*■ °°, which is of no 
import here except that it required the ruling out of solutions 
which grow in space. Because two "columns" of information are 
needed at the boundary to start the solution for g, another 
condition is required. Recalling the discussion in the previ- 
ous chapter concerning the modelling of boundary conditions, 
this problem of requiring an extra condition could lead one to 
expend great effort in attempting to develop a consistent and 
meaningful extra condition, when in fact if one simply applies 
the method of parametric differentiation as it is defined, the 
second condition will become apparent. The point is that the 
boundary is a real thing, not some abstraction in parameter 
space, and consequently the source excitation has no depen- 
dence on ip whatsoever; therefore one can in fact, without wor- 
rying about what £ is at the boundary (which, incidentally, 
would overspecify the problem) , write another boundary condi- 
tion. 
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g (0, t) = 0 (5.19) 

This completes the information needed to march out the solu- 
tion of the problem. It is interesting to note that because g 
and g are zero at the boundary, the only non-zero term in the 

X 

finite difference cube on its first application is the inhomo- 
geneous term. Without this term, there would be no solution 
because all terms in the finite difference cube would be zero. 
Again, this is consistent with intuitive thoughts about the 
form of the solution for g; it should in fact be very small 
near the boundary so that when the quadrature is performed to 
calculate the particle displacement, the net change in the 
value of £ from the inviscid to the viscous solution is small. 
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VI . RESULTS 

The practical realization of the preceding development is 
a computer program which solves the non-linear plane wave pro- 
pagation problem in a viscous fluid and provides output which 
allows one to readily compare the results obtained with this 
method, and those obtained by Keck and Beyer [4] . Their solu- 
tion was chosen as the point for comparison because it was 
simple and easy to calculate, and because this analysis had to 
be conducted on a shoestring budget. Unfortunately, this did 
not allow good comparisons for the higher ranges of T because 
the Keck/Beyer solution is not sufficiently accurate [6], 
lacking harmonic content beyond the sixth harmonic. One can 
argue though that the results for large T appear to be consis- 
tent with one's expectations, and that the method appears to 
generate a plausible progression of results as the value of T 
is decreased (corresponding to an increase in \]j) . An annotated 
version of this computer program appears in Appendix A. It is 
worthwhile to note at this point that because of the hyperbolic 
nature of the problem, the domain of dependence for a given 
point becomes larger and larger the further away from the 
boundary that the point is located. The result is that if one 
desires a solution value N discrete steps away from the boun- 
dary, N 2 + 2N + 1 points are required in the solution grid. 

This is a fundamental limitation of this method as it is pres- 
ently formulated, and may be difficult to improve upon. 
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In this analysis the Mach number was taken to be rather 
large, e = O(.Ol), in order to minimize N while still being 
able to generate useful results. B/A was taken as 6, which is 
the approximate value of that quantity for water; thus 3=4. 
For a given 3e product, the discontinuity distance, X, will 
occur at a point which is 1/2tt3e wavelengths from the boundary, 
regardless of frequency. Thus, for 3e = .04, X 4A. In water, 
a Mach number of . 01 corresponds to an acoustic pressure given 
by the expression p = p o C^e or p ^ 225 atmospheres, or in deci- 
bel units, SPL = 267 dB re 1 uPa, which is quite a strong sig- 
nal. The values of have been varied from zero to .01 which 
would reflect a change from inviscid to highly viscous propa- 
gation. Recalling the form ip = vb/C Q A , it is clear that by 
fixing ip , the frequency for which the solution is valid is al- 
so fixed for a given fluid. To put these results in a proper 
perspective, a value of = . 01 in water would correspond to a 
frequency of approximately five gigahertz. However, the re- 
sults for smaller values of ip represent more reasonable fre- 
quencies. Another interpretation of the results for ip = .01 
is that a solution has been obtained for the fluid with an ex- 
tremely high viscosity coefficient such as glycerin, a solu- 
tion which has not previously been obtained analytically, be- 
cause of the restriction a/K <<1, which is necessary to re- 
duce the equations to analytically tractable forms. One might 
conclude from these results that there is nothing significantly 
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different in the form of the solution for such a liquid. In 
any case, the results which follow will indicate the potential 
of this method for application over a great range of values of 
the non-linear and viscous parameters. 

Figures 5a, b, c, and d are plots of the particle dis- 
placement calculated by parametric differentiation, the par- 
ticle velocity calculated from this solution, the difference 
between the inviscid (Fubini) and the viscous (Keck-Beyer) 
displacement solutions, and the difference between the invis- 
cid and viscous (parametric differentiation - P.D.) displace- 
ment solutions, respectively, for e = .01, and \j) = .01 (F = 
1.27). The two plots of the calculated difference are the ba- 
sic results used for evaluation of the parametric differentia- 
tion solution. The differences are calculated according to 
Equation 4.16, and one would expect the calculated differences 
to be in phase with the displacement solutions; since viscos- 
ity should reduce the absolute amplitudes, the difference 
should be maximum when the Fubini solution has a maximum and 
minimum when the Fubini solution has a minimum. The zeroes 
will not necessarily coincide since the effect of the non-lin- 
earity is to lengthen the compressive portion of the waveform 
while reducing its magnitude and to shorten the rarefactive 
portion while increasing its absolute magnitude; viscosity will 
tend to keep the waveform sinusoidal, thus maintaining the 
length of the two portions of the waveform. Thus one could 
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hardly expect the zeroes to coincide. 

Since T = 1.27, one can infer that the Keck-Beyer solu- 
tion should be accurate, and also that some harmonic content 
should be evident in the particle velocity solution, assuming 
T = 1 is an accurate representation of the incipience of im- 
portant non-linear effects. In analyzing Figures 5c and 5d , 
one is immediately struck by the fact that although the calcu- 
lated differences seem to be of the same form, there seems to 
be some sort of ramp type behavior entering into the P.D. dis- 
placement solution. If one were to try to place a physical 
interpretation on this result, it would be that the ability of 
the medium to restore itself elastically has been exceeded due 
to the large displacement amplitudes and that therefore a per- 
manent displacement has been induced, just as if a spring had 
been stretched beyond its elastic limit inducing a permanent 
strain. Due to the fluid nature of the medium and because no 
analytic solution predicts such behavior, such a physical in- 
terpretation may not have much meaning. An attempt will be 
made to offer a possible explanation for this result in a la- 
ter paragraph, but this remains an unsolved problem in this 
work. The average numerical loss in the positive portion of 
the waveform is clearly greater than that of the negative por- 
tion, rather than being equal as one would expect. 

The velocity solution, on the other hand, seems to be 
quite normal. One can observe that though non-linear effects 
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are producing some steepening of the profile, viscous dissipa- 
tion does in fact reduce the amplitudes properly; the ampli- 
tude of the last peak is 3.1 dB down vice 3.3 dB down by the 
Keck-Beyer solution, a difference of about 6%. The last point 
in the velocity profile is numerically inaccurate because of 
the lack of sufficient points at the tip of the solution mesh. 
The small "bump" in the velocity profile is not easily ex- 
plained, but probably arises from numerical inaccuracies. It 
seems remarkable that with little direct effort aimed at re- 
ducing computational inaccuracies in the numerical scheme, the 
difference in the solutions is but 6%, remembering all the 
while that this difference may not represent an error at all, 
but rather a more accurate solution produced by a method which 
does not resort to common approximations. One may wonder at 
this point why an equation in particle displacement vice one 
in particle velocity was treated, the answer being that the 
equation in particle velocity would still have required one to 
know the displacement solution in order to arrive at the vari- 
able coefficients which would have appeared in the parameter 
space equation, and thus would probably have been no essential 
simplification. Another aspect of the P.D. solution which is 
worthy of note is the steepness of the velocity profile. The 
Keck-Beyer solution for this value of T looks essentially like 
a decaying sinusoid, while the P.D. solution appears to indi- 
cate that the non-linearity remains fairly important, result- 



- 60 - 



ing in some steepening of the profile. This suggests that 
perhaps T = 1 may not be a sufficiently strong indicator of 
the incipience of non-linear effects. 

Results analogous to those of Figure 5 have been obtained 
for different values of the problem parameters, but with the 
number of points utilized in the calculation maintained at 
N = .01. Figures 6 and 7 are the plotted results for e = .01 
and ^ = .001 (T = 12.7), and for e = .01 and \p = .0000463 (T 
= 275) respectively. Figures 8, 9, and 10 are the results for 
e = .02 and ip = .01 (T = 2.55), ip = .001 (T = 25.5), and ip = 
.0000463 (T = 550) respectively. In addition, calculations 
have been conducted for e = .0075 and e = .005, but are not 
shown here. These results are not much different in form from 
those discussed previously, but they uncover some problems with 
the solution technique as well as help demonstrate that the 
method is in fact a potentially useful one for investigation 
of non-linear wave propagation problems, since it successfully 
solves the equation of motion in a predictable and consistent 
manner for various values of T. The range of the 3e product 
had to be restricted as it was because of the need to minimize 
computer costs and thus the number of points at which the dif- 
ference equation was to be solved. 

Shifting attention to Figure 6, first recall that with T 
= 12.7 the Keck-Beyer solution is probably only marginally ac- 
curate. The displacement solution again exhibits the ramp type 
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behavior, although it is only readily clear when analyzing 
Figure 6d , which shows a greater average loss in the positive 
portion of the waveform than in the negative portion. Figure 
6c shows that the Keck-Beyer solution is in fact no longer ac- 
curate near the discontinuity distance. The velocity waveform 
appears to be consistent with expectations with only slight 
attenuation because of the relative smallness of the attenua- 
tion coefficient over the spatial length of the calculations. 
However, small signal attenuation, i.e., that attenuation which 
would occur if no harmonics were non-linearly produced, would 
predict a local Mach number of .0092 at the discontinuity while 
in this case a value of .0096 was obtained. This result is 
not consistent; one would expect greater attenuation in the 
non-linear case because the generation of harmonics will re- 
duce the magnitude of the primary in addition to small signal 
attentuation , while the harmonics themselves will be attenu- 
ated even more rapidly because of greater attenuation coeffi- 
cients. Unfortunately, in this case there is no easy explana- 
tion of the result; however, it did prompt some thoughts about 
what the Fubini solution looks like and what the changing 
shape of the waveform should imply with respect to the instan- 
taneous values of particle velocity near the discontinuity 
distance . 

From the point of view of energy conservation, one would 
expect that in a fluid without viscous dissipation, the time 
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or space average energy density in a plane wave train propagating 
isentropically will remain constant. The implication of this 
is that as the finite amplitude waveform steepens and eventu- 
ally assumes a sawtooth shape, the peak values of the particle 
velocity must increase in order to maintain constant average 
energy density in the wave train. For a sawtooth wave to have 
the same energy density as a sinusoid of unit amplitude and 
identical period, its amplitude must be approximately 1.22. 
Realizing that a fully formed sawtooth develops at a distance 
from the boundary of about 3 . 5X [3] , it would still be expected 
that the peak value of the particle velocity at the discontin- 
uity distance would be greater than that at the boundary. The 
base solution (Fubini) calculated in this study is accurate to 
five places in displacement and three in velocity, and the 
peak particle velocity at the discontinuity distance is ident- 
ical to that at the boundary, implying that there in fact has 
been an energy loss in an inviscid solution, and it is certain- 
ly not clear where this energy has gone. In any case, for a 
T of 12.7, non-linear effects are definitely important, and it 
is entirely conceivable that the P.D. solution has somehow 
recognized an error in the base solution and has yielded a 
quadrature value which correctly accounts for the steepening 
of the profile and thus predicts a particle velocity which is 
greater than that which would be obtained with a sinusoid and 
small signal attenuation. This explanation is of course sup- 
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position since there are too few points calculated to obtain 
a good space average energy density to verify an energy loss 
which is greater than that predicted by small signal attenua- 
tion. It indicates an unsolved problem and at the same time 
suggests that parametric differentiation predicts the proper 
solution to the problem. Blackstock [9] suggests the mechan- 
ism for dissipation in an inviscid fluid after a shock has 
formed; here, there is no reason for such dissipation since 
there are no shocks in the first zone of non-linear propaga- 
tion. 

Figure 7 shows much the same results as those of Figures 
5 and 6. In this case, T = 275, and the Keck-Beyer solution 
is clearly not converging to the proper result as evidenced by 
Figure 7c. Again the ramp behavior is clear in Figure 7d, 
though the magnitudes of the difference values are understand- 
ably smaller. The peak value of the calculated Mach number at 
the discontinuity distance is now .0106 which is in fact great- 
er than the peak at the boundary. The differences between the 
inviscid and viscous displacement solutions are now so small 
that on the scale of these plots, the viscous solution appears 
identical to the inviscid solution. This is not remarkable, 
but it is noteworthy that the plotted differences, though 
small, produce a smooth curve rather than one that jumps from 
point to point. Figures 8, 9, and 10, which represent the 
calculations for a source Mach number of .02, show results 
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which are similar to those shown in Figures 5, 6, and 7, and 
do not need to be analyzed with as much detail. However, in 
Figure 8, with r = 2.55, the peak Mach value near the discon- 
tinuity distance is .0147 vice .0137 according to the Keck- 
Beyer solution for a difference of 7%, which is consistent 
with the 6% difference obtained in the case e = .01, T = 1.27. 

It is also consistent with the idea that peak Mach numbers 
should increase in the non-linear case since small signal at- 
tenuation would predict a Mach number of .014 at the same 
point for a purely sinusoidal wave. The ramp behavior is evi- 
dent in each of the Figures 8d, 9d, and lOd. Finally, these results 
also provide a counter example for the idea that the peak 
Mach number should increase near the discontinuity distance in 
the case of strong non-linearity, i.e. when F is greater than 
10. The particle velocity plots of Figures 9b and 10b indi- 
cate peak Mach values which remain less than the peak at the 
boundary; in Figure 9b the value is .0193 and in Figure 10b it 
is .01998. Small signal attenuation would predict a value of 
.0192 for the result of Figure 9b. Energy considerations dic- 
tate an increase in the peak Mach values and in this case it 
does not occur. However, the spatial step size used to arrive 
at the plots of Figures 6 and 7 was .04, while in the calcula- 
tions used to obtain the plots of Figures 9 and 10, the spa- 
tial step size was .02. Recalling the finite difference ap- 
proximations utilized, this difference in step size may be the 
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explanation for the difference in results. The truncation er- 
ror for the second-order difference approximations is of O^Ax) 2 .^ 
Therefore, it is certainly conceivable that the step size of 
.04 is not accurate enough to yield good results for particle 
velocity. At the same time the considerations of energy can- 
not be ignored. It is likely that the increase in the peak 
Mach number in the results of Figures 6 and 7 is a result of 
numerical inaccuracy, but this is certainly not entirely clear. 
The other calculations conducted with Mach numbers of .0075 
and . 005 produced results which were consistent with the other 
results obtained, but they are not presented here because the 
spatial step size was sufficiently large that the results are 
somewhat questionable because of truncation error. 

There remains a question with respect to the curious ramp 
behavior in the particle displacement solutions. In examining 
the plots of the difference between the inviscid and viscous 
solutions, this behavior is always clear, but would not neces- 
sarily be clear at all in analyzing the particle displacement 
results, particularly as ^ becomes smaller and smaller. ip is 
a direct measure of the viscosity of the assumed fluid rather 
than a relative measure of viscosity such as T. Recalling that 
a f of .01 would imply a highly viscous liquid or a very high 
frequency wave and that the approximation upon which previous 
authors' results are based is a/K << 1, which implies weak 
viscosity, the reason for the ramp may become clear. Previous 
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analyses, in making this approximation, may have neglected a 
rather interesting physical interpretation of the non-linear 
loss mechanism. For small ip , the values of the calculated 
differences are so small compared to the values of particle 
displacement that this effect would not be noted when solv- 
ing for displacement or velocity without consideration of the 
differences between the viscous and inviscid solutions. In 
fact, in making the assumption a/K << 1, the values of ip are 
apparently restricted to relatively small values; a/K = 



1TV ^ > ; and ^ = f~t ' therefore it is clear that previous 
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analyses have not considered values of ip on the order of those 
considered here. For large ip and thus highly viscous propaga- 
tion, previous results are not valid, and therefore if the 
ramp behavior has a physical explanation, it could not have 
been noted before. 

In the parameter space Equations 4.6 or 4.9, there are 
four terms which involve the dependent variable / g , g , 
g , and g , which result from similar terms in the displace- 

XX u X 

ment space Equation 3.9. g is a term which relates to the 

XX 

elastic properties of the medium, g to the inertial proper- 
ties, and 9 xxt to the viscous properties. The g^ term implies 
that there exists some modification of the elastic properties 
of the medium. Recalling the discussion in the previous sec- 
tion concerning the effect of this term in the portions of the 
waveform where the velocity gradient u assumes large values, 

X 



- 67 - 



the source of the ramp behavior may become evident. When E, 

XX 

is large, the g term will dominate the equation. In the in- 

X 

viscid solution, E becomes infinite at the first shock. 

^xx 

Since the inviscid solution is in fact che starting solution 
and the quadrature from parameter space to displacement space 
is accomplished with an integration with respect to i|/, it is 
obviously not unreasonable to expect E, to assume large val- 
ues when ip is small; also, one would expect E, to assume nu- 
merical values which are more accurate for the current \ p step 
of the calculation if the spatial step size was small in the 
regions of large E, . In any case, if the g term dominates 

XX X 

the equation in a specific region, the g solution in that re- 
gion will vary linearly with x, hence the E, solution will vary 
in a similar manner, and hence the ramp. 

In the absence of viscosity, it is clear that the fluid 
particles oscillate about an equilibrium position; with vis- 
cosity, realizing that the g^ term is one which, in essence, 
alters the elastic constant of the system, it is postulated 
that the following phenomenon) occurs. Because of the large 
velocity amplitudes, the phase velocity increases in the com- 
pressive portion of the wave and decreases in the rarefactive 
portion; this is well known. These large particle velocities 
produce correspondingly large particle displacements, but now 
the fluid particles, instead of being purely elastic, do not 
have all of their kinetic energy at zero displacement converted 
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to potential energy at full displacement. Some is lost to 
small signal attenuation of both fundamental and harmonics, 
and the rest produces an irreversible displacement of the 
fluid particles; that is, work is expended on the fluid system 
to produce a displacement of the equilibrium positions of the 
fluid particles. Unfortunately, this effect can only be pos- 
tulated here because a precise energy balance has not been 
performed. In order to calculate the specific energy lost in 
producing a given displacement at a given point, the net ac- 
celeration on the fluid particles which produced this displace- 
ment would be required, and this calculation has not yet been 
performed . 

The reason for this apparent displacement of the equili- 
brium positions of the fluid particles may be explained in a 
more rhetorical manner as an interaction between non-linear 
and viscous effects. The non-linear effect is to alter the 
elastic constant thus changing the speed of propagation of the 
wave, resulting in the steepening of the waveform. Now, be- 
cause the oncoming compressive portion of a given cycle of the 
wave is travelling faster than the preceding rarefactive por- 
tion and because viscosity is acting to slow the velocity of 
the fluid particles in the rarefactive portion even more, the 
particles, after passage of the given cycle, do not return to 
their former equilibrium position, but are instead minutely 
displaced. The oncoming compressive portion of the wave is 
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' filling' the preceding rarefactive portion and consequently 
the fluid particles in the rarefactive portion do not have the 
opportunity to recover their initial equilibrium positions. 

Such an effect would obviously tend to create a vacuum behind 
the wave train, but this can be justified by requiring that 
fluid flow in from outside the region excited by the source. 

Finally, there are some matters concerning the numerical 
method which require attention. The stability criterion de- 
rived in the preceding chapter was arrived at in a rather un- 
conventional manner. Some numerical experimentation was per- 
formed to determine just how sensitive to the ratio of step 
sizes the solution was. The exact stability criterion was 

Ax < At/(1 + If) 26 (6.1) 

which implies for 9£/3x positive, c = .01, and 3 = 4 that Ax 
_< At/1.08 near the first shock. This should be an absolute 
limit on the spatial step size since 3£/3x will have its maxi- 
mum value at this point. The equations were solved for Ax=.04 
and At=.04, and the numerical solution became highly unstable 
beyond the second cycle of the waveform. Then At was set equal 
to .045 which would be just slightly greater than the At re- 
quired by the stability criterion for Ax=.04, and the solution 
was stable throughout. All subsequent calculations were con- 
ducted with a convenient At selected so that it was large enough 
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for the given Ax to satisfy the stability criterion. 

Another matter of concern in the calculations is the ef- 
fect of the step sizes, in x and t and in t|j as well, on the 
accuracy of the solution. In this analysis, the base solution 
was calculated accurate to five significant figures, obviously 
limiting the accuracy of the calculated solution to the same 
extent. The base solution required 51 harmonics at the dis- 
continuity distance for this accuracy. The spatial step sizes 
used for the various solutions were Ax = .08, .06, .04, and 

.02. Since the truncation errors in the second order differ- 
ences are of O^Ax) 2 ), it would be reasonable to assume that the 
calculated values of g would be accurate to two or three places 
depending on the step size used. However, to return to dis- 
placement space, the values of g are multiplied by a difference 
between two values of ip , and hence for a maximum ip of .01, if 
g is accurate to three significant figures, the differences 
calculated in the quadrature will leave the displacement solu- 
tion accurate to five significant figures, neglecting errors 
introduced by the discretization of the solutions with respect 
to ip . This then would be consistent with a base solution cal- 
culated to four significant figures, for Ax=.02 and marginally 
consistent for Ax=.04. This is the reason why the results for 
e = .005 and e = .0075 were not presented; their accuracy could 
only be considered marginal at best for the maximum ^ value. 
These were the solutions which employed step sizes of .06 and 
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.08. Roundoff error was not considered to be of significant 
concern in this analysis with machine accuracy to seven signi- 
ficant figures. The solutions also did not seem to be sensi- 
tive to the step size in time as long as the time step was 
large enough to satisfy the stability criterion at all points. 
Finally, the step size in ip seemed to matter very little in 
the form of the results if not in the exact numerical results. 
It is obvious that in integrating across orders of magnitude, 
the integration must be handled carefully, but experience dem- 
onstrated that values of g at a given point changed little as 
ip was varied, and thus the error introduced by utilizing a 
relatively large step size in ip should be small. Consequently, 
the following scheme was adopted: 

1) the first 'p value was the limiting value 

2) the second ip value was that which when used in 
the quadrature yielded difference values which 
were just within machine roundoff accuracy 

3) each subsequent step in the order of magnitude 
of \p was broken into an algebraically increasing 
number of steps on a logarithmic basis. 

For example, with i p Q = 0, \p^ = 10 ^ , there would be two steps 

between ip = 10 ^ and ip = 10 , with \p ^ = .316 x 10 ^ and ip^ 

-5 -5 -4 

= 10 , three steps between ^ = 10 and ^ = 10 , ip^ = .214 

-4 -4 -4 

x 10 , \p c = .463 x 10 , and if>, = 10 , etc. There is no 

5 D 



- 72 - 



clear method to determine the proper step size in i p, but this 
one produced consistent results. In general, it would be in- 
appropriate if g varied in a highly non-linear manner with ij >. 

The final comments concerning the numerical results are 
made with respect to the formulation of the problem. Equations 
4.6 and 4.9 were the parameter space equations developed from 
the equation treated by Keck and Beyer and from the parent e- 
quation without approximation, respectively. Early calcula- 
tions were performed with Equation 4.6 and later ones with E- 
quation 4.9. For e = .01 and = .01, the maximum percentage 
difference between the calculated solutions by Equations 4.6 
and 4.9 was approximately 4%, demonstrating that the equation 
solved in previous work does in fact produce very little error 
for Mach numbers up to .01. Finally, the two forms of the 
boundary values for g appeared to make little difference in 
the final result; producing a maximum percentage difference in 
the calculated values of particle displacement of approximate- 
ly 1% between the case utilizing the 'natural' boundary condi- 
tions for g and the case in which artificial dependence on ip 
was introduced at the boundary. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 

A theoretical study of non-linear plane acoustic wave 
propagation in a viscous fluid has been conducted utilizing a 
little known mathematical technique, parametric differentiation, 
in order to demonstrate its utility when applied to acoustic 
propagation problems, and to obtain whatever physical insight 
into the problem the results may present. It has been demon- 
strated that the method can in fact yield results which are 
consistent with other analytic solutions to the problem for 
certain arbitrarily selected values of the viscous and non- 
linear parameters. No restriction was required in the selec- 
tion of these values by the mathematical nature of the method, 
and there is no reason to believe that the method could not be 
applied for any values of the parameters. Restrictions which 
may exist in the application of the method are computational 
ones, specifically: 1) the results can only be as accurate as 

the solution to the problem for ip - ^ o ; 2) accuracy may be 
further limited by the finite difference approximations util- 
ized; 3) accuracy may be limited by the ability of the computer 
to solve the finite difference scheme for a large number of 
points with a given time and memory budget; for instance, dou- 
ble precision will require twice the computer memory; and 4) 
the hyperbolic nature of the problem requires calculations to 
be conducted throughout the domain of dependence of a given 
point to arrive at a solution at that point. These restric- 
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tions may render the method unusable for calculations involv- 
ing large propagation distances, but they do not detract from 
its ability to provide physical insight. Further research may 
provide a reasonable technique for application of the method 
when large distances are involved. The solution of the prob- 
lem by parametric differentiation has allowed the elimination 
of the restriction a/k << 1 heretofore assumed in typical pro- 
pagation studies. Elimination of this restriction has produced 
results which are valid for highly viscous fluids as well as 
demonstrated that the results for such fluids are essentially 
similar to those for less viscous fluids. It has allowed the 
formulation of a postulated behavior for all viscous fluids, 
which may explain physically why there are non-linearly in- 
duced energy losses above and beyond those associated with 
viscous attenuation. Specifically the solutions produced by 
this method predict that some of the energy contained in a 
finite amplitude wave in a viscous medium is expended in pro- 
ducing a displacement of the equilibrium positions of the fluid 
particles. These results have also been obtained cheaply; for 
101 spatial points, the computer time required to produce a 
result for single values of ^ and the Be product is approxi- 
mately .15 minutes with 400 kilobytes of computer memory. 

This is not to suggest that this work proves beyond a 
doubt that the method is a useful one. Consideration of prop- 
agation once a shock has actually formed has not been attempted. 
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If the shocks were in fact discontinuities as predicted by the 
inviscid theory, the technique probably would not be useful 
beyond the first zone of propagation. However, since the wave- 
form is in fact continuous through a shock of finite thickness, 
it is reasonable to expect that with denser sampling of the 
waveform through the shocks, good numerical results could be 
obtained. The next step in a continuing study of the solution 
of the problem through parametric differentiation would be to 
select Blackstock's weak shock solution [3], which connects 
the Fubini solution [1] with the Fay solution [2] in its in- 
viscid limit to yield a solution valid to the end of the sec- 
ond zone of propagation, as the base solution for ip Q = 0, and 
to attempt to extend the solution into the region of periodic 
shocks. With this result a precise energy balance calculation 
would demonstrate whether or not energy is lost in producing a 
permanent displacement of the fluid particles. 

More fundamentally, it would be useful if a mathematical 
analysis of the technique were to be conducted in order to de- 
termine why it apparently correctly predicts solutions for many 
different classes of problems, and what if any are the theor- 
etical limitations to its application. Beyond this, there are 
more interesting problems in non-linear w ave propagation to 
which the method could be applied. The first, and one which 
can be readily adapted to the computer program already written, 
would be an investigation of the solution of the Burgers equa- 
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tion as r is varied. Having noted previously that this solu- 
tion is difficult to calculate as T becomes much greater than 
x/X or for T > 50 because of slow convergence, it is suggested 
that the analytic solution be calculated for T = 1. Then with 
the Burgers equation in the form solved by Blackstock [5] , 
with V = u/U , 

V a - VVy = r_lv yy ' v (°' y) = sinwy (7.1) 

differentiate with respect to r to obtain, 

9 0 ■ Vg y - V = r " lg yy ' r " 2v yy’ g(o ' y) = 0 (7 ' 2) 

where a is a spacelike, and y a timelike variable. This equa- 
tion would be perfectly suited to solution on the mesh used in 
this analysis because there would be no value at the i+l,j-l 
point, and the scheme would be completely explicit. The solu- 
tion for F = 1 could then be extended to arbitrary values of T with- 
out significant computational problems. Because of the space 
scale of Equation 7.1, this technique will allow easy calcula- 
tion of the solution in both the first and second zones of 
propagation. 

Finally, there are the two problems of finite amplitude 
wave propagation which are yet to be solved analytically, those 
of propagation of cylindrical and spherical waves in viscous 
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fluids; parametric differentiation may provide an invaluable 
tool in solving these problems. The fundamental equations of 
motion do not lend themselves to analytic solution, there be- 
ing no simple transform such as that used by Blackstock in the 
plane wave problem to reduce the equation of motion to a Bur- 
gers equation [18] . Indeed, the inviscid forms of the equa- 
tion have no exact analytic solutions which have been published. 
Blackstock [9] has obtained approximate solutions to the cylin- 
drical and spherical wave problems for an inviscid fluid, one- 
dimensional wave propagation, and r >> 1, which are very sim- 
ilar in form to the Fubini solution; these results were ob- 
tained through manipulations of the equation of motion which 
yielded an inviscid Burgers equation. Such solutions could be 
used as one base solution for a solution to the viscous prob- 
lem by parametric differentiation, realizing that they are on- 
ly approximate. Other authors [12] [19] have produced particu- 
lar solutions for spherical and cylindrical wave propagation, 
but with an accurate base solution, parametric differentiation 
could be used to generate many solutions, for a fixed Be pro- 
duct, with relative economy. Consequently, the apparent course 
of study to pursue for non-linear cylindrical and spherical 
wave propagation is the development of improved base solutions 
for = 0 and the application of parametric differentiation 
to these solutions, producing a large number of solutions which 
may then allow for some generalization of the results. Such a 
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research program could lead to a significant gain in the under 
standing of the physical mechanisms which are collectively des 



cribed as non-linear attenuation in viscous fluids. 
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Figure 3 The Nine Point Finite Difference Cube 




Figure 4 The Solution Mesh, the Arrows Indicate 
the Marching Direction Utilized Herein 
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Figure 5a. Displacement Profile, Be = .04, r = 1.27 
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Figure 5b. Velocity Profile, = .04, T = 1.27 
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Figure 5d. Difference Between Fubini and P.D. Displacement Profiles 
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Figure 6a. Displacement Profile, 6e = .04, r = 12. 
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Figure 6b. Velocity Profile, 3e = .04, F = 12. 
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Figure 6c. Difference Between Fubini and Keck-Beyer Displacement Profiles 

= .04, T = 12.7 
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Figure 6d. Difference Between Fubini and P.D. Displacement Profiles 

8s = .04, r = 12.7 
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Figure 7a. Displacement Profile, Be = .04, r = 275 
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Figure 7b. Velocity Profile, Be = .04, T = 275 
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Figure 7c - Difference Between Fubini and Keck-Beyer Displacement Profiles 

Be = .04, r = 275 



- 100 - 




l .onu3iy3wnN-iNi9nd=dii30 



Figure 7d. Difference Between Fubini and P.D. Displacement Profiles 

ge = .04, r = 275 
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Figure 8c. Difference Between Fubini and Keck-Beyer 
Displacement Profiles, = .08, I’ = 2.55 
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Figure 8d. Difference Between Fubini and P.D. 
Displacement Profiles, 8e = .08, F = 2.55 
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Figure 9a 
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Figure 9b. Velocity Profile, Be = .08, r = 25.5 
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9c. Difference Between Fubini and Keck-Beyer 
Displacement Profiles, Be = .08, T = 25.5 
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Figure 9d. Difference Between Fubini and P.D. 
Displacement Profiles, Be = .08, F = 25.5 
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Figure 10b. Velocity Profile, $£ = .08, r = 550 
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Figure 10c. Difference Between 
Displacement Profiles, 3e 
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Figure lOd. Difference Between Fubini and P.D. 
Displacement Profiles, 8e = .08, r = 550 
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